function [err, tc_int]=TC_INT(sf, plot_flag)
% Analyze changes in intrinsic plasticity over time
% This returns a modified control worksheet
% Analysis results are placed in the sheet, and the data are all
% corrected in advance

% return structure is tc_int (time course of intrinsic changes), and contains:
% tc_int.timestamp = time of latest analysis
% tc_int.time[]  is the time point for each record
% tc_int.vm[] is the resting membrane potential for each record
% tc_int.fsl[] is the first spike latency for each record
% tc_int.fisi[] is the first interspike interval for each record
% tc_int.nspk[] is the spike count during s3
% tc_int.bridge[] is the estimated bridge balance
%



err = 0;
QueMessage('TC_INT - Starting');
dat = [];
time = [];
tc_int = []; % initialize it
do_spike = 1;	% 0 turns off detailed spike analysis; 1 turns it on

if(nargin == 0)
    sf = getmainselection;
    plot_flag = 1;
end;
if(nargin == 1)
    plot_flag = 1;
end;
%try

if(sf > 0)
    pflag = getplotflag;
    QueMessage('TC_INT analysis', 1); % clear the que
    for i = 1:length(sf)
        tc_int2(sf(i), pflag, do_spike);
    end;
end;



function tc_int2(sf, pflag, do_spike);
global VOLTAGE CURRENT DFILE
global CONTROL
% make sure data and selection are synchronized
[DFILE, DP, err] = ana_setup(DFILE, sf);

if(err == 1)
    fprintf(1, 'TC_INT: failed to open %s\n', CONTROL(sf).filename);
    return; % just return
end;
plot_flag = 1;
switch(DP.protocol)
    case 'drugmonitor'
        % abstract general information
        % mean times:
        m1=DP.tsw2+DP.tstepdur*0.6;
        m2=DP.tsw2+DP.tstepdur*0.99;
        m3=DP.tsw2+DP.tstepdur*0.97; % use for h current assessment

        %peak times:
        p1=DP.tsw2+DP.tstepdur*0.01;
        p2=DP.tsw2+DP.tstepdur*0.50;
        p3=DP.tsw2+DP.tstepdur*0.25;

    otherwise
        % abstract general information
        % mean times:
        m1=DP.tdel+DP.tstepdur*0.6;
        m2=DP.tdel+DP.tstepdur*0.99;
        m3=DP.tdel+DP.tstepdur*0.97; % use for h current assessment

        %peak times:
        p1=DP.tdel+DP.tstepdur*0.01;
        p2=DP.tdel+DP.tstepdur*0.50;
        p3=DP.tdel+DP.tstepdur*0.25;
end
% calculate measurement times

trmp = floor(0.95*DP.tdel./DP.RATES); % for RMP/ihold determination

mean_t1=floor(m1./DP.RATES); % for ss iv
mean_t2=floor(m2./DP.RATES);
mean_t3=floor(m3./DP.RATES);


peak_t1=floor(p1./DP.RATES); % for peak iv
peak_t2=floor(p2./DP.RATES);
peak_tfit=floor(p3./DP.RATES);

sw1 = floor(DP.tsw2./DP.RATES);
sw2 = floor(DP.tsw3./DP.RATES);
ssw1 = DP.tsw2;
ssw2 = DP.tsw3;
% compute the time base for plotting
tmax=max(max(DP.time));
% DFILE.ztime = find(DFILE.ztime > 0);

tc_int(sf).time = DP.TM;

% measure holding current and "rmp"
for rec = 1:DP.records
    tc_int(sf).ihold(rec) = mean(CURRENT(rec,1:trmp(rec))); % save in array for later usage
    tc_int(sf).vm(rec) = mean(VOLTAGE(rec,1:trmp(rec)));  % mean rmp each record
    % measure steady-state IV and voltage
    tc_int(sf).vpk(rec) = mean(VOLTAGE(rec,peak_t1(rec):peak_t2(rec)));
    tc_int(sf).vss(rec) = mean(VOLTAGE(rec,mean_t1(rec):mean_t2(rec)));
    tc_int(sf).istep(rec) = mean(CURRENT(rec,mean_t1(rec):mean_t2(rec))) - tc_int(sf).ihold(rec);
end
CONTROL(sf).Rmp = mean(tc_int(sf).vm); % mean across all records...
CONTROL(sf).iHold = mean(tc_int(sf).ihold); % across all records
tc_int(sf).rin = 1000 * (tc_int(sf).vss - tc_int(sf).vm) ./ tc_int(sf).istep; % compute input resistance over time

% do a time constant fit on the trace at -175 to -225 pA injection...
%QueMessage('TC_INT - Fitting membrane TC');
iH = CONTROL(sf).iHold;
a0 = -60;
a1 = -20;
tau = 10;
%[fpar, chisq, niter, volt_fit1(i,:)] = mrqfit('exponential', [a0 a1 tau], time_fit, vsmo(tfit(i), t_list), [], ...
%   [], [], [], 50, []);

bl = find(DP.ZT <= DP.cond_baseline);
% get spike times
QueMessage('TC_INT - Finding spikes');
[first_spike, first_isi, nr_spikes, spike_train]=find_spikes(DFILE, ssw1, ssw2, DP.spike_thresh);
if(~isempty([spike_train]))
    % disp(sprintf('Rec: found %4d spikes', sum(nr_spikes)))
    tc_int(sf).fisi = first_isi;
    tc_int(sf).fsl = first_spike;
    for rec=1:DP.records
        index = find([spike_train.source] == rec); % find out which one goes with this rec
        if(isempty(index))
            tc_int(sf).nspike(rec) = NaN; % no spikes, no index!
        else
            tc_int(sf).nspike(rec) = length(spike_train(index).latency);
        end;
    end;
    [fslp] = msd(first_spike, bl);
    [fisip] = msd(first_isi, bl);
else
    tc_int(sf).fisi = NaN*ones(1,DP.records);
    tc_int(sf).fsl = NaN*ones(1,DP.records);
    tc_int(sf).nspike = zeros(1,DP.records);

    [fslp] = [];
    [fisip] = [];
    
end;

[vmp] = msd(tc_int(sf).vm, bl);
tx=[min(DP.ZT),max(DP.ZT)];


CONTROL(sf).tc_int = tc_int; % make it a substructure - that's how we store results
QueMessage('TC_INT - Analysis complete');

% plot if figure is set
if plot_flag>0
    t0 = tc_int(sf).time(1);
    fsize = 7;
    msize = 3;
    tmax = max(DP.time);
    msg_pos = [0.37 0.00 0.15 0.07];

    h = findobj('Tag', 'TC_INT'); % check for pre-existing window
    if(isempty(h)) % if none, make one
        h = figure('Tag', 'TC_INT', ...
            'MenuBar', 'none', ...
            'WindowButtonMotionFcn', 'datac(''mouse_motion'', gcbf);', ...
            'WindowButtonDownFcn', 'datac(''mouse_down'', gcbf);', ...
            'WindowButtonUpFcn', 'datac(''mouse_up'', gcbf);', ...
            'NumberTitle', 'off');
    end
    figure(h); % otherwise, select it
    clf; % always clear the window...
    %Command Menu
    uimenu('Label', 'Close &Window', 'Position', 1, 'callback', 'close(findobj(''Tag'', ''TC_INT''));' );
    uimenu('Label', '&Print', 'Callback', 'print;');
    uimenu('Label', 'Print&Setup', 'Callback', 'printdlg;');
    set(h, 'Name', sprintf('Intrinsic Time Course Analysis - File: %s', CONTROL(sf).filename))

    subplot('Position', [0.1, 0.90, 0.8, 0.07]);
    axis([0,1,0,1]);
    axis('off');

    text(0.35,0.7,sprintf('%s  R: [%d-%d]',CONTROL(sf).filename, ...
        number_arg(CONTROL(sf).recbeg), number_arg(CONTROL(sf).recend)), 'FontSize', fsize);
    text(0.40,0.3,sprintf('%s', date),'FontSize', fsize);

    subplot('Position', [0.1 0.50 0.8 0.4]);
    plot(tc_int(sf).time-t0, tc_int(sf).fsl, 'rx', 'MarkerSize', 5, 'MarkerFaceColor', 'r')
    hold on
    plot(tc_int(sf).time-t0, tc_int(sf).fisi, 'bo', 'MarkerSize', 4, 'MarkerFaceColor', 'b')
    if(~isempty(fisip))
        pmsd(tx, fisip, 'b');
    end;

    if(~isempty(fslp))
        pmsd(tx, fslp, 'r');
    end;
    hold off
    set(gca, 'FontSize', fsize, 'XTickLabel', '');
    xlabel('');
    ylabel('fsl, fisi');
    datac_setcrosshair(gca, 'TC_INT_info', 'mn', 'ms', msg_pos);

    subplot('Position', [0.1 0.4 0.8 0.08]);
    plot(tc_int(sf).time-t0, tc_int(sf).nspike, 'ko', 'MarkerSize', 4, 'MarkerFaceColor', 'k');
    [yl] = ylim;
    yln(1)=0; yln(2)=max(tc_int(sf).nspike)+1;
    ylim(yln);
    set(gca, 'FontSize', fsize, 'XTickLabel', '');
    xlabel('');
    ylabel('# spikes');
    datac_setcrosshair(gca, 'TC_INT_info', 'mn', 'N', msg_pos);

    subplot('Position', [0.1 0.3 0.8 0.08]);
    plot(tc_int(sf).time-t0, tc_int(sf).ihold, 'ko', 'MarkerSize', 4, 'MarkerFaceColor', 'k');
    [yl] = ylim;
    if(abs(yl(1)-yl(2))<50)
        yln(1)=mean(yl)-25; yln(2)=mean(yl)+25;
        ylim(yln);
    end;
    xlabel('');
    set(gca, 'FontSize', fsize, 'XTickLabel', '');
    ylabel('I_h');
    datac_setcrosshair(gca, 'TC_INT_info', 'mn', 'Ih', msg_pos);

    subplot('Position', [0.1 0.2 0.8 0.08]);
    plot(tc_int(sf).time-t0, tc_int(sf).rin, 'ko', 'MarkerSize', 4, 'MarkerFaceColor', 'k');
    [yl] = ylim;
    if(abs(yl(1)-yl(2))<10)
        yln(1)=mean(yl)-5; yln(2)=mean(yl)+5;
        ylim(yln);
    end;
    xlabel('');
    set(gca, 'FontSize', fsize, 'XTickLabel', '');
    ylabel('R_in');
    datac_setcrosshair(gca, 'TC_INT_info', 'mn', 'Mohm', msg_pos);

    subplot('Position', [0.1 0.1 0.8 0.08]);
    plot(tc_int(sf).time-t0, tc_int(sf).vm, 'ko', 'MarkerSize', 4, 'MarkerFaceColor', 'k');
    [yl] = ylim;
    if(abs(yl(1)-yl(2))<20)
        yln(1)=mean(yl)-10; yln(2)=mean(yl)+10;
        ylim(yln);
    end;
    hold on;
    pmsd(tx, vmp, 'k');
    hold off
    set(gca, 'FontSize', fsize);
    ylabel('Vm');
    datac_setcrosshair(gca, 'TC_INT_info', 'mn', 'mV', msg_pos);

    orient landscape

    % control printing and closing of window for automatic runs
    % plot_flag = 1 creates plot and leaves it up
    % plot_flag = 2 creates plot but closes it when done
    % plot_flag = 2 creates plot and prints it and then closes it
    if (plot_flag == 2)
        print -dljet3;
    end
    if plot_flag >= 2
        close
    end
end

function [p] = msd(x, l)
p = [];
bl = find(~isnan(x(l))); % only for ones with valid data...
if(~isempty(bl))
    m = mean(x(bl));
    s = std(x(bl));
    p.m = [m,m];
    p.p = [p.m+s]; % for stddev lines
    p.n = [p.m-s];
end;
return;

function pmsd(tx, p, c)
plot(tx, p.m, 'LineStyle', '-', 'color', c);
plot(tx, p.p, 'LineStyle', '--', 'color', c);
plot(tx, p.n, 'LineStyle', '--', 'color', c);
return;

